Abundancia Relativa

Instalar y cargar paqueterías

if(!requireNamespace("tidyverse", quietly=TRUE))
  install.packages("tidyverse")
library(tidyverse)
library(plotly)

Importar los datos a analizar

De las tablas de abundancia de los distintos amplicones, generadas con el procesamiento de dada2, se importan los archivos generados.

data <- read.delim("~/Downloads/secuencias/Tablas/ASVs_total_silva.csv") #Adecuar la ruta del archivo

Obtener abundancia relativa de NA

Para tener una idea del porcentaje de lecturas que van a ser modificadas debido a la falta de categorización taxonómica, se obtiene la abundancia de NAs a lo largo de todas las muestras y ASVs.

# Sumar el conteo de NA y el total de lecturas para todas las muestras
total_NA <- data %>%
  filter(is.na(Family)) %>%
  select(starts_with("AC")) %>%
  pivot_longer(everything(), names_to = "Sample", values_to = "Count") %>%
  summarize(Total_NA = sum(Count))

# Sumar el total de lecturas en todas las muestras
total_reads <- data %>%
  select(starts_with("AC")) %>%
  pivot_longer(everything(), names_to = "Sample", values_to = "Count") %>%
  summarize(Total_Reads = sum(Count))

# Calcular la abundancia relativa de NA en todo el conjunto de muestras
abundancia_relativa_NA <- total_NA$Total_NA / total_reads$Total_Reads * 100

# Mostrar el resultado
abundancia_relativa_NA

Calcular abundancia relativa

Familia

Para calcular la abundancia relativa, primero se deben reemplazar los valores nulos NAs por la categoría “No Clasificados”. Posteriormente se calcula la abundancia relativa por familia taxonómica.

data$Family <- ifelse(is.na(data$Family), "No Clasificados", as.character(data$Family))

family_abundance <- data %>%
  select(starts_with("AC"), Family) %>%
  pivot_longer(cols = -Family, names_to = "Sample", values_to = "Count") %>%
  group_by(Sample, Family) %>%
  summarize(Total = sum(Count)) %>%
  mutate(Percentage = Total / sum(Total) * 100)

A los ASVs con una abundancia menor al 1% se les asigna en la nueva categoría de Otros.

family_abundance <- family_abundance %>%
  mutate(Family = ifelse(Percentage < 1, "Otros", Family))

Posteriormente se procede a realizar de nuevo el cálculo de abundancia relativa.

family_abundance <- family_abundance %>%
  group_by(Sample, Family) %>%
  summarize(Total = sum(Total)) %>%
  mutate(Percentage = Total / sum(Total) * 100)

Género

Para los taxones de Género y Especie se sigue un paso adicional. Se reemplazan los valores nulos NA con el nombre del taxón Familia usando el formato “u_Nombre”.

##Se usa de nuevo la tabla original sin la categoría de Otros para no perder esas asignaciones taxonómicas 
data <- read.delim("~/Downloads/secuencias/Tablas/ASVs_total_silva.csv") 
## Reemplazar los valores NA en la columna "Genus" usando el taxon anterior.
data$Genus <- ifelse(is.na(data$Genus), paste0("u_", data$Family), data$Genus)
## Reemplazar "u_NA" y "u_Otros" con NA en la columna "Genus"
data$Genus <- ifelse(data$Genus == "u_NA" | data$Genus == "u_Otros", NA, data$Genus)
# Reemplazar los valores NA en la columna "Genus" con "No Clasificados"
data <- data %>%
  mutate(Genus = if_else(is.na(Genus), "No Clasificados", Genus))
genus_abundance <- data %>%
  select(starts_with("AC"), Genus) %>%
  pivot_longer(cols = -Genus, names_to = "Sample", values_to = "Count") %>%
  group_by(Sample, Genus) %>%
  summarize(Total = sum(Count)) %>%
  mutate(Percentage = Total / sum(Total) * 100)
genus_abundance <- genus_abundance %>%
  mutate(Genus = ifelse(Percentage < 1, "Otros", Genus))
genus_abundance <- genus_abundance %>%
  group_by(Sample, Genus) %>%
  summarize(Total = sum(Total)) %>%
  mutate(Percentage = Total / sum(Total) * 100)

Especies

Por último se puede realizar lo mismo con el taxón de especies para mayor especificidad.

##Se usa de nuevo la tabla original sin la categoría de Otros para no perder esas asignaciones taxonómicas 
data <- read.delim("~/Downloads/secuencias/Tablas/ASVs_total_silva.csv") 
## Reemplazar los valores NA en la columna "Species" usando el taxon de familia.
data$Species <- ifelse(is.na(data$Species), paste0("u_", data$Family), data$Species)
## Reemplazar "u_NA" y "u_Otros" con NA en la columna "Species"
data$Species <- ifelse(data$Species == "u_NA" | data$Species == "u_Otros", NA, data$Species)
# Reemplazar los valores NA en la columna "Species" con "No Clasificados"
data <- data %>%
  mutate(Species = if_else(is.na(Species), "No Clasificados", Species))
species_abundance <- data %>%
  select(starts_with("AC"), Species) %>%
  pivot_longer(cols = -Species, names_to = "Sample", values_to = "Count") %>%
  group_by(Sample, Species) %>%
  summarize(Total = sum(Count)) %>%
  mutate(Percentage = Total / sum(Total) * 100)
species_abundance <- species_abundance %>%
  mutate(Species = ifelse(Percentage < 1, "Otros", Species))
species_abundance <- species_abundance %>%
  group_by(Sample, Species) %>%
  summarize(Total = sum(Total)) %>%
  mutate(Percentage = Total / sum(Total) * 100)

Graficar abundancias

Para poder visualizar las abundancias relativas, se procede a graficar los diferentes taxones. En las gráficas generadas, Otros encapsula los ASVs con taxonomía no encontrada en la base de entrenamiento, y aquellos cuya abundancia relativa era menor al 1%.

Familia

family_abundance$Family <- factor(family_abundance$Family, 
                                  levels = c(setdiff(unique(family_abundance$Family), "Otros"), "Otros"))
family_plot<-ggplot(family_abundance, aes(x = Sample, y = Percentage, fill = Family)) +
  geom_bar(stat = "identity", position = "stack") +
  labs(title = "Abundancia Relativa por Familia", y = "Porcentaje", x = "Muestra") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1), 
        legend.key.size = unit(0.5, "cm"),
        legend.text = element_text(size = 8),
        plot.margin = margin(1, 1, 1, 1, "cm")
        ) + scale_color_brewer(palette="Set3") + guides(fill = guide_legend(ncol = 2))
ggplotly(family_plot)

Género

genus_abundance$Genus <- factor(genus_abundance$Genus, 
                                  levels = c(setdiff(unique(genus_abundance$Genus), "Otros"), "Otros"))
genus_plot<-ggplot(genus_abundance, aes(x = Sample, y = Percentage, fill = Genus)) +
  geom_bar(stat = "identity", position = "stack") +
  labs(title = "Abundancia Relativa por Género", y = "Porcentaje", x = "Muestra") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1),
       legend.key.size = unit(0.5, "cm"),
        legend.text = element_text(size = 8),
        plot.margin = margin(1, 1, 1, 1, "cm")
       ) +scale_color_brewer(palette="Set3") + guides(fill = guide_legend(ncol = 2))
ggplotly(genus_plot)

Especies

species_abundance$Species <- factor(species_abundance$Species, 
                                  levels = c(setdiff(unique(species_abundance$Species), "Otros"), "Otros"))
species_plot<-ggplot(species_abundance, aes(x = Sample, y = Percentage, fill = Species)) +
  geom_bar(stat = "identity", position = "stack") +
  labs(title = "Abundancia Relativa por Especies", y = "Porcentaje", x = "Muestra") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1),
        legend.key.size = unit(0.5, "cm"),
        legend.text = element_text(size = 8),
        plot.margin = margin(1, 1, 1, 1, "cm")
       ) + scale_color_brewer(palette="Set3") + guides(fill = guide_legend(ncol = 2))

ggplotly(species_plot)

Guardar gráficas

Para guardar las gráficas en distintos formatos, se tiene que asignar una variable a la función ggplot. A continuación se muestra un ejemplo de como quedarían las líneas de código para realizar esto.

path<-"/Users/adriantorres/Downloads/secuencias"
ggsave(file.path(path, "family_abundance.svg"), plot = family_plot, dpi = 600)
ggsave(file.path(path,"family_abundance.jpg"), plot = family_plot, dpi = 600)
ggsave(file.path(path,"family_abundance.pdf"), plot = family_plot, dpi = 600)
#En este caso family_plot corresponde al nombre de la variable asignada a la función ggplot. 

Se realiza lo mismo con el taxón género y especies.

ggsave(file.path(path, "genus_abundance.svg"), plot = genus_plot, dpi = 600)
ggsave(file.path(path,"genus_abundance.jpg"), plot = genus_plot, dpi = 600)
ggsave(file.path(path,"genus_abundance.pdf"), plot = genus_plot, dpi = 600)
ggsave(file.path(path, "species_abundance.svg"), plot=species_plot, dpi = 600)
ggsave(file.path(path,"species_abundance.jpg"), plot=species_plot, dpi = 600)
ggsave(file.path(path,"species_abundance.pdf"), plot=species_plot, dpi = 600)